home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_01 / bnarth2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  6.9 KB  |  277 lines
  1. /*
  2.  *   MIRACL arithmetic routines 2 - multiplying and dividing BIG NUMBERS.
  3.  *   bnarth2.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include "miracl.h"
  8. #define abs(x)  ((x)<0? (-(x)) : (x))
  9. #define sign(x) ((x)<0? (-1) : 1)
  10.  
  11. /* Access global variables */
  12.  
  13. extern small base;    /* number base     */
  14. extern int  nib;      /* length of bigs  */
  15. extern int  depth;    /* error tracing ..*/
  16. extern int  trace[];  /* .. mechanism    */
  17. extern big  w0;       /* workspace big   */
  18. extern bool check;    /* error checking  */
  19.  
  20. void multiply(x,y,z)
  21. big x;
  22. big y;
  23. big z;
  24. {  /*  multiply two big numbers: z=x.y  */
  25.     int i,sz,xl,yl,j,ti;
  26.     small carry;
  27.     if (ERNUM) return;
  28.     if (size(y)==0 || size(x)==0) 
  29.     {
  30.         zero(z);
  31.         return;
  32.     }
  33.     depth++;
  34.     trace[depth]=5;
  35.     if (TRACER) track();
  36.     if (notint(x) || notint(y))
  37.     {
  38.         berror(12);
  39.         depth--;
  40.         return;
  41.     }
  42.     sz=sign(x[0])*sign(y[0]);
  43.     xl=abs(x[0]);
  44.     yl=abs(y[0]);
  45.     if (yl==1)
  46.     {
  47.         premult(x,y[1],z);
  48.         insign(sz,z);
  49.         depth--;
  50.         return;
  51.     }
  52.     zero(w0);
  53.     if (xl+yl>nib && (check || (xl+yl)>2*nib))
  54.     {
  55.         berror(3);
  56.         depth--;
  57.         return;
  58.     }
  59.     if (x==y)
  60.     { /* squaring can be done nearly twice as fast */
  61.         for (i=1;i<xl;i++)
  62.         { /* long multiplication */
  63.             carry=0;
  64.             for (j=i+1;j<=xl;j++)
  65.             { /* Only do above the diagonal */
  66.                 carry=muldiv(x[i],x[j],w0[i+j-1]+carry,base,&w0[i+j-1]);
  67.             }
  68.             w0[xl+i]=carry;
  69.         }
  70.         w0[0]=xl+xl-1;
  71.         padd(w0,w0,w0);     /* double it */
  72.         carry=0;
  73.         for (i=1;i<=xl;i++)
  74.         { /* add in squared elements */
  75.             ti=i+i;
  76.             carry=muldiv(x[i],x[i],w0[ti-1]+carry,base,&w0[ti-1]);
  77.             w0[ti]+=carry;
  78.             carry=0;
  79.             if (w0[ti]>=base)
  80.             {
  81.                 carry=1;
  82.                 w0[ti]-=base;
  83.             }
  84.         }
  85.     }
  86.     else for (i=1;i<=xl;i++)
  87.     { /* long multiplication */
  88.         carry=0;
  89.         for (j=1;j<=yl;j++)
  90.         { /* multiply each digit of y by x[i] */
  91.             carry=muldiv(x[i],y[j],w0[i+j-1]+carry,base,&w0[i+j-1]);
  92.         }
  93.         w0[yl+i]=carry;
  94.     }
  95.     w0[0]=(xl+yl)*sz; /* set length and sign of result */
  96.     lzero(w0);
  97.     copy(w0,z);
  98.     depth--;
  99. }
  100.  
  101. void divide(x,y,z)
  102. big x;
  103. big y;
  104. big z;
  105. {  /*  divide two big numbers  z=x/y : x=x mod y  *
  106.     *  returns quotient only if  divide(x,y,x)    *
  107.     *  returns remainder only if divide(x,y,y)    */
  108.     small carry,borrow,try,r,ldy,ra,d,tst,dig,pdiff,psum;
  109.     int i,k,m,sz,sx,sy,x0,y0,w00;
  110.     if (ERNUM) return;
  111.     depth++;
  112.     trace[depth]=6;
  113.     if (TRACER) track();
  114.     if (x==y) berror(7);
  115.     if (notint(x) || notint(y)) berror(12);
  116.     if (size(y)==0) berror(2);
  117.     if (ERNUM)
  118.     {
  119.         depth--;
  120.         return;
  121.     }
  122.     sx=sign(x[0]);    /* extract signs ... */
  123.     sy=sign(y[0]);
  124.     sz=sx*sy;
  125.     x[0]=abs(x[0]);  /* ... and force operands to positive */
  126.     y[0]=abs(y[0]);
  127.     x0=x[0];
  128.     y0=y[0];
  129.     copy(x,w0);
  130.     w00=w0[0];
  131.     if (check && (w00-y0+1>nib))
  132.     {
  133.         berror(3);
  134.         depth--;
  135.         return;
  136.     }
  137.     d=0;
  138.     if (x0==y0)
  139.     {
  140.         if (x0==1) /* special case - x and y are both ints */
  141.         { 
  142.             d=w0[1]/y[1];
  143.             w0[1]%=y[1];
  144.             lzero(w0);
  145.         }
  146.         else if (w0[x0]/4<y[x0])
  147.         while (compare(w0,y)>=0)
  148.         {  /* small quotient - so do up to four subtracts instead */
  149.             psub(w0,y,w0);
  150.             d++;
  151.         }
  152.     }
  153.     if (compare(w0,y)<0)
  154.     {  /*  x less than y - so x becomes remainder */
  155.         if (x!=z)  /* testing parameters */
  156.         {
  157.             copy(w0,x);
  158.             insign(sx,x);
  159.         }
  160.         if (y!=z) convert(sz*d,z);
  161.         insign(sy,y);
  162.         depth--;
  163.         return;
  164.     }
  165.     if (y0==1)
  166.     {  /* y is small - so use subdiv instead */
  167.         r=subdiv(w0,y[1],w0);
  168.         if (y!=z)
  169.         {
  170.             copy(w0,z);
  171.             insign(sz,z);
  172.         }
  173.         if (x!=z)
  174.         {
  175.             convert(r,x);
  176.             insign(sx,x);
  177.         }
  178.         insign(sy,y);
  179.         depth--;
  180.         return;
  181.     }
  182.     if (y!=z) zero(z);
  183.     d=base/(y[y0]+1);   /* have to do it the hard way */
  184.     premult(w0,d,w0);
  185.     premult(y,d,y);
  186.     ldy=y[y0];
  187.     for (k=w00;k>=y0;k--)
  188.     {  /* long division */
  189.         if (w0[k+1]==ldy) /* guess next quotient digit */
  190.         {
  191.             try=base-1;
  192.             ra=ldy+w0[k];
  193.         }
  194.         else try=muldiv(w0[k+1],base,w0[k],ldy,&ra);
  195.         while (ra<base)
  196.         {
  197.             tst=muldiv(y[y0-1],try,(small)0,base,&r); 
  198.             if (tst< ra || (tst==ra && r<=w0[k-1])) break;
  199.             try--;  /* refine guess */
  200.             ra+=ldy;
  201.         }
  202.         m=k-y0;
  203.         if (try>0)
  204.         { /* do partial subtraction */
  205.             borrow=0;
  206.             for (i=1;i<=y0;i++)
  207.             {
  208.               borrow=muldiv(try,y[i],borrow,base,&dig);
  209.               pdiff=w0[m+i]-dig;
  210.               if (pdiff<0)
  211.               { /* set borrow */
  212.                   borrow++;
  213.                   pdiff+=base;
  214.               }
  215.               w0[m+i]=pdiff;
  216.             }
  217.             if (w0[k+1]<borrow)
  218.             {  /* whoops! - over did it */
  219.                 w0[k+1]=0;
  220.                 carry=0;
  221.                 for (i=1;i<=y0;i++)
  222.                 {  /* compensate for error ... */
  223.                     psum=w0[m+i]+y[i]+carry;
  224.                     carry=0;
  225.                     if (psum>=base)
  226.                     {
  227.                         carry=1;
  228.                         psum-=base;
  229.                     }
  230.                     w0[m+i]=psum;
  231.                 }
  232.                 try--;  /* ... and adjust guess */
  233.             }
  234.             else w0[k+1]-=borrow;
  235.         }
  236.         if (k==w00 && try==0) w00--;
  237.         else if (y!=z) z[m+1]=try;
  238.     }
  239.     if (y!=z) z[0]=(w00-y0+1)*sz; /* set sign and length of result */
  240.     w0[0]=y0;
  241.     lzero(y);
  242.     lzero(z);
  243.     if (x!=z)
  244.     {
  245.         lzero(w0);
  246.         subdiv(w0,d,x);
  247.         insign(sx,x);
  248.     }
  249.     subdiv(y,d,y);
  250.     insign(sy,y);
  251.     depth--;
  252. }
  253.  
  254. void mad(x,y,z,w,q,r)
  255. big x,y,z,w,q,r;
  256. { /* Multiply, Add and Divide; q=(x*y+z)/w remainder r   *
  257.    * returns remainder only if w=q, quotient only if q=r *
  258.    * add done only if x, y and z are distinct.           */
  259.     if (ERNUM) return;
  260.     depth++;
  261.     trace[depth]=24;
  262.     if (TRACER) track();
  263.     check=OFF;           /* turn off some error checks */
  264.     if (w==r)
  265.     {
  266.         berror(7);
  267.         depth--;
  268.         return;
  269.     }
  270.     multiply(x,y,w0);
  271.     if (x!=z && y!=z)add(w0,z,w0);
  272.     divide(w0,w,q);
  273.     if (q!=r) copy(w0,r);
  274.     check=ON;
  275.     depth--;
  276. }
  277.